process fgmax output

After running the GeoClaw code, use this notebook to read in and plot the fgmax results.

In [1]:
%matplotlib notebook
In [2]:
from pylab import *
In [3]:
from clawpack.geoclaw import fgmax_tools, topotools, dtopotools
from clawpack.visclaw import colormaps
from scipy.interpolate import RegularGridInterpolator
import numpy
In [4]:
GEmap = imread('EagleHarborGE.png')
GEextent = [-122.55,-122.48,47.61,47.64]
In [5]:
figure(figsize=(8,6))
imshow(GEmap,extent=GEextent)
gca().set_aspect(1./cos(47.6*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20);
In [6]:
fg = fgmax_tools.FGmaxGrid()
In [7]:
fg.read_input_data('fgmax_eagle_harbor.txt')
In [8]:
fg.read_output(outdir='_output')

Reading _output/fort.FG1.valuemax ...
Reading _output/fort.FG1.aux1 ...

Define colormaps:

In [9]:
zmin = -60.
zmax = 40.
land_cmap = colormaps.make_colormap({ 0.0:[0.1,0.4,0.0],
                                     0.25:[0.0,1.0,0.0],
                                      0.5:[0.8,1.0,0.5],
                                      1.0:[0.8,0.5,0.2]})

sea_cmap = colormaps.make_colormap({ 0.0:[0,0,1], 1.:[.8,.8,1]})

cmap, norm = colormaps.add_colormaps((land_cmap, sea_cmap),
                                     data_limits=(zmin,zmax),
                                     data_break=0.)

Compute subsidence during earthquake

The values stored in fg.B are the cell-averaged topography values used in GeoClaw after uplift or subsidence due to the earthquake. Since Eagle Harbor lies right on top of the Seattle Fault, there is significant subsidence in this region. It is important to correct for this before plotting the "maximum flow depth on shore", for example, since it affects what cells are considered "on shore", as we illustrate below. We generally want to show the flooding depth in regions that are on shore before the earthquake, even if they end up below sea level afterward.

Read in the sea floor deformation file and interpolate to the locations of the fgmax grid. Note this is not exactly the right correction in cases like this where the uplift varies greatly over a small area, but is pretty good and fine for most cases.

In [10]:
dtopo_path = '../dtopo/dtopofiles/seattlefault_uniform.tt3'
dtopo = dtopotools.DTopography(dtopo_path, dtopo_type=3)
x1d = dtopo.X[0,:]
y1d = dtopo.Y[:,0]
dtopo_func = RegularGridInterpolator((x1d,y1d), dtopo.dZ[-1,:,:].T, 
                method='linear', bounds_error=False, fill_value=0.)
        
dz = dtopo_func(list(zip(numpy.ravel(fg.X), numpy.ravel(fg.Y))))
dz = numpy.reshape(dz, fg.X.shape)

# Estimate of topography before earthquake:
B0 = fg.B - dz

print('The maximum subsidence over the fgmax region is %.2f meters' % dz.min())
print('The maximum uplift over the fgmax region is %.2f meters' % dz.max())

The maximum subsidence over the fgmax region is -1.78 meters
The maximum uplift over the fgmax region is 6.90 meters

There is massive uplift just south of Eagle Harbor but subsidence on the harbor shore, as seen in this plot of the deformation (with the original coast in green):

In [11]:
figure(figsize=(8,6))
ax = axes()
contour(fg.X, fg.Y, B0, [0], colors='g')  # original shoreline
dtopo.plot_dZ_colors(t=2, axes=ax)
axis(GEextent)
gca().set_aspect(1./cos(47.6*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20);

Compare topography before and after event

In [12]:
figure(figsize=(8,10))
subplot(2,1,1)
pc = pcolormesh(fg.X, fg.Y, B0, cmap=cmap, norm=norm)  
cb = colorbar(pc, extend='both', shrink=0.5)
cb.set_label('meters')
gca().set_aspect(1./cos(47.6*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20);
title('GeoClaw topography B, before earthquake')

subplot(2,1,2)
pc = pcolormesh(fg.X, fg.Y, fg.B, cmap=cmap, norm=norm)  
cb = colorbar(pc, extend='both', shrink=0.5)
cb.set_label('meters')
gca().set_aspect(1./cos(47.6*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20);
title('GeoClaw topography B, after earthquake')
Out[12]:
<matplotlib.text.Text at 0xb1dac5588>

Compare the shorelines

In [13]:
figure(figsize=(8,6))
imshow(GEmap,extent=GEextent)
contour(fg.X, fg.Y, B0, [0], colors='b')
contour(fg.X, fg.Y, fg.B, [0], colors='r')

gca().set_aspect(1./cos(47.6*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20);
title('Shore before (blue) and after (red) earthquake')
Out[13]:
<matplotlib.text.Text at 0xb1fc4d860>

Define onshore depth based on B0

Note that you get misleading results if you replace B0 by fg.B below!

In [14]:
h_onshore = ma.masked_where(B0 < 0., fg.h)
print('Maximum flow depth onshore: %.2f m' % h_onshore.max())

Maximum flow depth onshore: 12.94 m
In [15]:
#bounds_depth = array([1e-6,0.25,0.5,0.75,1,1.25,1.5])
bounds_depth = array([1e-6,0.5,1.,2.,3.,4.,5.])

cmap_depth = mpl.colors.ListedColormap([[.7,.7,1],\
                 [.5,.5,1],[0,0,1],\
                 [1,.7,.7], [1,.4,.4], [1,0,0]])

# Set color for value exceeding top of range to purple:
cmap_depth.set_over(color=[1,0,1])

# Set color for land points without inundation to light green:
cmap_depth.set_under(color=[.7,1,.7])

norm_depth = mpl.colors.BoundaryNorm(bounds_depth, cmap_depth.N)
    
figure(figsize=(8,6))
pc = pcolormesh(fg.X, fg.Y, h_onshore, cmap=cmap_depth, norm=norm_depth)
cb = colorbar(pc, extend='max', shrink=0.5)
cb.set_label('meters')
gca().set_aspect(1./cos(48*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20)
title('Maximum flow depth')
#savefigp('h_onshore.png')
Out[15]:
<matplotlib.text.Text at 0xb1e198e10>

Plot on map image

To plot on a map we want to mask out the onshore dry points so they don't all get colored green and obscure the map.

(Again, observe what happens is you replace B0 by fg.B).

In [16]:
h_onshore = ma.masked_where(B0 < 0, fg.h)

# define h_wet_onshore to suppress plotting dry onshore fgmax points:
h_wet_onshore = ma.masked_where(fg.h < 0.01, h_onshore)

figure(figsize=(8,6))
imshow(GEmap,extent=GEextent)

pc = pcolormesh(fg.X, fg.Y, h_wet_onshore, cmap=cmap_depth, norm=norm_depth)
cb = colorbar(pc, extend='max', shrink=0.5)
cb.set_label('meters')
gca().set_aspect(1./cos(48*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20)
title('Maximum flow depth')
Out[16]:
<matplotlib.text.Text at 0xb1ea81400>

Plot the maximum speed

In [17]:
print('Maximum flow speed: %.2f m/s' % fg.s.max())

Maximum flow speed: 11.90 m/s
In [18]:
bounds_speed = np.array([1e-6,0.5,1.5,2,2.5,3,4.5,6])
cmap_speed = mpl.colors.ListedColormap([[.9,.9,1],[.6,.6,1],\
                 [.3,.3,1],[0,0,1],\
                 [1,.7,.7], [1,.4,.4], [1,0,0]])

# Set color for value exceeding top of range to purple:
cmap_speed.set_over(color=[1,0,1])

# Set color for land points without inundation to light green:
cmap_speed.set_under(color=[.7,1,.7])

#imshow(GEmap,extent=GEextent)
norm_speed = mpl.colors.BoundaryNorm(bounds_speed, cmap_speed.N)

figure(figsize=(8,6))
pc = pcolormesh(fg.X, fg.Y, fg.s, cmap=cmap_speed, norm=norm_speed)
cb = colorbar(pc, extend='max')
cb.set_label('m/s')

contour(fg.X, fg.Y, B0, [0], colors='k')  # plot original shoreline

gca().set_aspect(1./cos(48*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20)
title('Maximum speed')
Out[18]:
<matplotlib.text.Text at 0xb1eb6fac8>

Plot on the map

Mask out the speed where the depth is 0 before plotting over the map.

We also include a transparency alpha in the call to pcolormesh so you can see the image below.

In [19]:
s_wet = ma.masked_where(fg.h < 0.01, fg.s)

figure(figsize=(8,6))
imshow(GEmap,extent=GEextent)

pc = pcolormesh(fg.X, fg.Y, s_wet, alpha = 0.2, 
                cmap=cmap_speed, norm=norm_speed)
cb = colorbar(pc, extend='max')
cb.set_label('m/s')
gca().set_aspect(1./cos(48*pi/180.))
ticklabel_format(useOffset=False)
xticks(rotation=20)
title('Maximum speed')
Out[19]:
<matplotlib.text.Text at 0xb1ee4c550>
In [20]:
close('all')